Skip to main content
Open in GitHub
Run in Google Colab

Introduction to Bluebook for Bulldozers Dataset

Overview of the Dataset​

The Bluebook for Bulldozers dataset, sourced from Kaggle, is a comprehensive compilation of historical auction prices for used bulldozers. It includes detailed information about the equipment, such as model type, usage hours, and age, providing an in-depth view of the market dynamics for these heavy machines.

Importance of the Dataset​

This dataset is particularly valuable for understanding price trends and market behavior in the heavy machinery sector. It serves as a rich source for predictive modeling, allowing analysts to forecast future prices and trends in the construction equipment market.

Explainable Machine Learning

Understanding Explainability in ML​

Explainable Machine Learning (ML) refers to techniques and methods that make the output of ML models more understandable to humans. This is crucial in scenarios where decision-making needs to be transparent and justifiable.

Role of Explainability in Predictive Modeling​

In predictive modeling, especially with complex models like those used in the Bluebook dataset, explainability helps in understanding how and why certain predictions are made. This transparency is vital for trust and validation of model predictions in real-world applications.

Potential Use Case: Analysing Arbitrage Opportunities​

Defining Arbitrage in Equipment Sales​

Arbitrage, in the context of equipment sales, refers to the opportunity to profit from price discrepancies in different markets or time periods. By accurately predicting the future prices of bulldozers, one can identify undervalued machines to purchase and overvalued machines to sell, capitalising on these market inefficiencies.

Utilizing ML Predictions for Arbitrage​

This subsection explores how the predictions made by the explainable ML model can be used to identify potential arbitrage opportunities in the bulldozer market. The explainable aspect of the model allows for a deeper understanding of the factors influencing price predictions, aiding in more strategic buying and selling decisions.

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_absolute_percentage_error, r2_score, mean_squared_error, explained_variance_score, mean_squared_log_error
import matplotlib.pyplot as plt
import xplainable as xp
from xplainable.core.ml.regression import XRegressor
from xplainable.core.optimisation.genetic import XEvolutionaryNetwork
from xplainable.core.optimisation.layers import Evolve, Tighten
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import warnings
import os
# from xplainable.visualisation.regression import plot_error
print(f"This notebook was created using Xplainable version {xp.__version__}")
Out: 
This notebook was created using Xplainable version 1.2.3

Error plot Function Overwrite​

Note: This will be part of the python package as part of v1.1 patch

from sklearn.metrics import mean_absolute_error
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

def plot_error(model, x, y, alpha=0.5, color_column=None):
    fig, ax = plt.subplots(figsize=(12, 8))

    y_pred = model.predict(x)
    mae = mean_absolute_error(y, y_pred)
    errors = abs(y - y_pred)

    if color_column is not None:
        if color_column not in x.columns:
            raise ValueError(f"The color_column {color_column} is not in the DataFrame.")
        
        # Convert column to categorical and get codes and unique values
        categories = x[color_column].astype('category').cat.categories
        codes = x[color_column].astype('category').cat.codes
        unique_codes = np.unique(codes)
        scatter = ax.scatter(y, y_pred, c=codes, alpha=alpha, cmap='plasma')

        # Create a legend with the actual category labels
        handles = [plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=scatter.cmap(scatter.norm(code)), 
                               markersize=10) for code in unique_codes]
        ax.legend(handles, categories, title=color_column)

    else:
        scatter = ax.scatter(y, y_pred, c=errors, alpha=alpha, cmap='plasma')
        plt.colorbar(scatter, ax=ax, label='Absolute Error')

    # Line for perfect predictions
    max_val = np.maximum(y.max(), y_pred.max())
    ax.plot([0, max_val], [0, max_val], 'k--', lw=2)

    # Labels and title
    ax.set_xlabel('True Values')
    ax.set_ylabel('Predicted Values')
    ax.set_title(f'Scatter Plot of True vs Predicted Values with MAE: {mae:.2f}')

    plt.show()

# Example usage:
# plot_error(model, X_train, y_train, alpha=0.5, color_column='ModelID')

Building the Dataset

It's possible to download the Bluebook dozer price prediction dataset at the following link: https://www.kaggle.com/c/bluebook-for-bulldozers/data

Following extraction of the .zip file build the dataset as below:

df = pd.read_csv('data/TrainAndValid.csv', parse_dates = ['saledate'])
df.head()
SalesIDSalePriceMachineIDModelIDdatasourceauctioneerIDYearMadeMachineHoursCurrentMeterUsageBandsaledate...Undercarriage_Pad_WidthStick_LengthThumbPattern_ChangerGrouser_TypeBackhoe_MountingBlade_TypeTravel_ControlsDifferential_TypeSteering_Controls
011392466600099908931571213200468Low2006-11-16...nannannannannannannannanStandardConventional
111392485700011765777121319964640Low2004-03-26...nannannannannannannannanStandardConventional
21139249100004348087009121320012838High2004-02-26...nannannannannannannannannannan
31139251385001026470332121320013486High2011-05-19...nannannannannannannannannannan
411392531100010573731731112132007722Medium2009-07-23...nannannannannannannannannannan

Add the machine appendix to concatenate information about the dozer assets

ma = pd.read_csv('data/Machine_Appendix.csv')
ma.head()
MachineIDModelIDfiModelDescfiBaseModelfiSecondaryDescfiModelSeriesfiModelDescriptorfiProductClassDescProductGroupProductGroupDescMfgYearfiManufacturerIDfiManufacturerDescPrimarySizeBasisPrimaryLowerPrimaryUpper
01131355350L350nannanLHydraulic Excavator, Track - 50.0 to 66.0 Metr...TEXTrack Excavators199426CaterpillarWeight - Metric Tons5066
14343538416C416CnannanBackhoe Loader - 14.0 to 15.0 Ft Standard Digg...BLBackhoe Loaders199726CaterpillarStandard Digging Depth - Ft1415
25343538416C416CnannanBackhoe Loader - 14.0 to 15.0 Ft Standard Digg...BLBackhoe Loaders199826CaterpillarStandard Digging Depth - Ft1415
37183538416C416CnannanBackhoe Loader - 14.0 to 15.0 Ft Standard Digg...BLBackhoe Loaders200026CaterpillarStandard Digging Depth - Ft1415
417531580D5GLGPD5GnanLGPTrack Type Tractor, Dozer - 85.0 to 105.0 Hors...TTTTrack Type Tractors200626CaterpillarHorsepower85105

Merging the dataset on the MachineID to extract useful information:

  • Find the columns that exist within the machine dictionary that aren't in the training dataset
  • Merge the new columns on the existing train dataset to enrich the information
new_cols = [col for col in ma.columns if col not in df.columns]
ma[new_cols].head()
MfgYearfiManufacturerIDfiManufacturerDescPrimarySizeBasisPrimaryLowerPrimaryUpper
0199426CaterpillarWeight - Metric Tons5066
1199726CaterpillarStandard Digging Depth - Ft1415
2199826CaterpillarStandard Digging Depth - Ft1415
3200026CaterpillarStandard Digging Depth - Ft1415
4200626CaterpillarHorsepower85105
merge_col = "MachineID"
df = pd.merge(df, ma[[merge_col]+ new_cols ],on='MachineID', how='left')

Feature Engineering Overview

Note: The approach to feature engineering presented here is foundational and does not encompass the full depth typically seen in data science projects. The release of version 1.1 has improved our system's ability to handle missing values directly, which reduces the time needed to process and analyse data that contains null entries.

# Extract Purchase date information
df['saleyear'] = df['saledate'].dt.year
df['salemonth'] = df['saledate'].dt.month
df['saledayofweek'] = df['saledate'].dt.day_name()

# Drop the Sale Date following extraction of features
df.drop('saledate', inplace=True, axis=1)
#Rename columns
df.columns = [col.replace("_","") for col in df.columns]
#Filter our erroneous purchase values
df = df[df.YearMade > 1920]
#Turn model id into a categorical value so it doesn't create regression splits
df["ModelID"] = df["ModelID"].astype(str) 
df
SalesIDSalePriceMachineIDModelIDdatasourceauctioneerIDYearMadeMachineHoursCurrentMeterUsageBandfiModelDesc...SteeringControlsMfgYearfiManufacturerIDfiManufacturerDescPrimarySizeBasisPrimaryLowerPrimaryUppersaleyearsalemonthsaledayofweek
0113924666000.099908931571213.0200468.0Low521D...Conventional2004.025CaseHorsepower110.0120.0200611Thursday
1113924857000.0117657771213.019964640.0Low950FII...Conventional1996.026CaterpillarHorsepower150.0175.020043Friday
2113924910000.043480870091213.020012838.0High226...nan2001.026CaterpillarOperating Capacity - Lbs1351.01601.020042Thursday
3113925138500.010264703321213.020013486.0HighPC120-6E...nan2010.0103KomatsuHorsepower225.0250.020115Thursday
4113925311000.01057373173111213.02007722.0MediumS175...nan2007.0121BobcatOperating Capacity - Lbs1601.01751.020097Thursday
..................................................................
412693633334410000.01919201214351492.02005nannan30NX...nan2005.02552IHIWeight - Metric Tons2.03.020123Wednesday
412694633334510500.01882122214361492.02005nannan30NX2...nan2005.02552IHIWeight - Metric Tons3.04.020121Saturday
412695633334712500.01944213214351492.02005nannan30NX...nan2005.02552IHIWeight - Metric Tons2.03.020121Saturday
412696633334810000.01794518214351492.02006nannan30NX...nan2006.02552IHIWeight - Metric Tons2.03.020123Wednesday
412697633334913000.01944743214361492.02006nannan30NX2...nan2005.02552IHIWeight - Metric Tons2.03.020121Saturday

Train on the top 6 dozers assets by count​

For timeliness of training filter the data on the Top 6 assets by count

models_to_train = df.ModelID.value_counts().index[:10].to_list()
models_to_train
Out: 
['4605',
 '3538',
 '4604',
 '3170',
 '3362',
 '3537',
 '4603',
 '3171',
 '3357',
 '3178']
df[df.ModelID.isin(models_to_train)]
SalesIDSalePriceMachineIDModelIDdatasourceauctioneerIDYearMadeMachineHoursCurrentMeterUsageBandfiModelDesc...SteeringControlsMfgYearfiManufacturerIDfiManufacturerDescPrimarySizeBasisPrimaryLowerPrimaryUppersaleyearsalemonthsaledayofweek
5113925526500.0100127446051213.02004508.0Low310G...nan2004.043John DeereStandard Digging Depth - Ft14.015.0200812Thursday
10113927824000.0102499846051213.020041414.0Medium310G...nan2004.043John DeereStandard Digging Depth - Ft14.015.020088Thursday
15113929119000.0100481046041213.019992450.0Medium310E...nan1999.043John DeereStandard Digging Depth - Ft14.015.0200611Thursday
62113946923000.0105886931711213.019989987.0High580L...nan1998.025CaseStandard Digging Depth - Ft14.015.020075Thursday
82113951533000.0101556546051213.020021268.0Medium310G...nan2002.043John DeereStandard Digging Depth - Ft14.015.020047Thursday
..................................................................
410243628823918200.01835461460414999.0200048.0Low310E...nan2000.043John DeereStandard Digging Depth - Ft14.015.020122Wednesday
410244628824025250.0190391446051490.020051988.0Low310G...nan2005.043John DeereStandard Digging Depth - Ft14.015.020121Saturday
410245628824125250.01860549460514999.02006nannan310G...nan2006.043John DeereStandard Digging Depth - Ft14.015.020124Wednesday
410246628824325000.0184618446051491.02006nannan310G...nan2006.043John DeereStandard Digging Depth - Ft14.015.020123Thursday
410264628834620500.0186708746041494.02000nannan310E...nan2000.043John DeereStandard Digging Depth - Ft14.015.020122Monday
data = df[df.ModelID.isin(models_to_train)]
m = data.isna().sum()
data = data[[col for col in data.columns if col not in data.columns[m == len(data)]]]

#Drop cols cardinality of 1
s = data.nunique()
car_cols = data.columns[(s == 1)]
data = data.drop(columns=car_cols)

#Update numeric columns to be float64
n_cols = data.select_dtypes(include=np.number).columns.tolist()
data[n_cols] = data[n_cols].astype('float64')
data.head()
SalesIDSalePriceMachineIDModelIDdatasourceauctioneerIDYearMadeMachineHoursCurrentMeterUsageBandfiModelDesc...TireSizeMfgYearfiManufacturerIDfiManufacturerDescPrimarySizeBasisPrimaryLowerPrimaryUppersaleyearsalemonthsaledayofweek
51.13926e+06265001.00127e+06460512132004508Low310G...nan200443John DeereStandard Digging Depth - Ft1415200812Thursday
101.13928e+06240001.025e+064605121320041414Medium310G...nan200443John DeereStandard Digging Depth - Ft141520088Thursday
151.13929e+06190001.00481e+064604121319992450Medium310E...nan199943John DeereStandard Digging Depth - Ft1415200611Thursday
621.13947e+06230001.05887e+063171121319989987High580L...nan199825CaseStandard Digging Depth - Ft141520075Thursday
821.13952e+06330001.01556e+064605121320021268Medium310G...nan200243John DeereStandard Digging Depth - Ft141520047Thursday

Addressing Multicollinearity in Model Interpretability​

It's well-understood in data science that multicollinearity can significantly hamper the interpretability of models, particularly those based on linear assumptions. The code snippet above demonstrates a rudimentary approach to mitigating multicollinearity by removing highly correlated features. However, it's important to acknowledge that this is a simplified illustration; in practice, the interplay between features can be more subtle and complex.

For robust feature selection and to enhance model explainability, we employ automated feature selection techniques that are thoroughly documented in our project's documentation. These methods go beyond pairwise correlations, considering the multidimensional structure of the data to retain the most informative features. While the current example is not exhaustive, it serves to highlight a fundamental step in preprocessing for linear models. Practitioners are encouraged to leverage our automatic feature selection capabilities to refine their models further and to ensure that the explanatory variables employed are truly reflective of independent factors influencing the response variable.

data["AgeAtSale"] = df["saleyear"] - df["MfgYear"]
drop_cols = [
    # "saleyear", #--> Data encoded in Age at Sale
    "MfgYear", #--> Data encoded in Age at Sale
    "YearMade", #--> Multicollinearity with MfgYear
            ]
target = 'SalePrice'
id_columns=["SalesID",'MachineID','auctioneerID','datasource']

Split the train and validation set​

data_train = data[data.saleyear!=2012]
data_val = data[data.saleyear==2012]
data_train= data_train.drop(columns=drop_cols)
data_val=data_val.drop(columns=drop_cols)

#Create the training and validation set
X_train, y_train = data_train.drop('SalePrice', axis=1), data_train['SalePrice']
X_valid, y_valid = data_val.drop('SalePrice', axis=1), data_val['SalePrice']
model = XRegressor(ignore_nan=False)

Initial fit of the Regressor​

model.fit(X_train, y_train,id_columns=id_columns)
Out: 
<xplainable.core.ml.regression.XRegressor at 0x284eff040>
model.evaluate(X_train, y_train)
Out: 
{'Explained Variance': 0.8476,
 'MAE': 4292.3689,
 'MAPE': 0.1616,
 'MSE': 39079631.4865,
 'RMSE': 6251.3704,
 'RMSLE': nan,
 'R2 Score': 0.8476}
model.explain()

Optimising the Model​

network = XEvolutionaryNetwork(model)

# Add the layers
# Start with an initial Tighten layer
network.add_layer(
    Tighten(
        iterations=100,
        learning_rate=0.1,
        early_stopping=20
        )
    )

# Add an Evolve layer with a high severity
network.add_layer(
    Evolve(
        mutations=100,
        generations=50,
        max_severity=0.5,
        max_leaves=20,
        early_stopping=20
        )
    )

# Add another Evolve layer with a lower severity and reach
network.add_layer(
    Evolve(
        mutations=100,
        generations=50,
        max_severity=0.3,
        max_leaves=15,
        early_stopping=20
        )
    )

# Add a final Tighten layer with a low learning rate
network.add_layer(
    Tighten(
        iterations=100,
        learning_rate=0.025,
        early_stopping=20
        )
    )

# Fit the network (before or after adding layers)
network.fit(X_train.drop(columns=id_columns), y_train)

# Run the network
network.optimise()
Out: 
  0%|          | 0/100 [00:00<?, ?it/s]
  0%|          | 0/50 [00:00<?, ?it/s]
  0%|          | 0/50 [00:00<?, ?it/s]
  0%|          | 0/100 [00:00<?, ?it/s]
<xplainable.core.optimisation.genetic.XEvolutionaryNetwork at 0x2d32d0730>
model.evaluate(X_train, y_train)
Out: 
{'Explained Variance': 0.8446,
 'MAE': 4134.8991,
 'MAPE': 0.1486,
 'MSE': 40127897.0119,
 'RMSE': 6334.6584,
 'RMSLE': nan,
 'R2 Score': 0.8435}

Simply by fitting a combination of 6 Tighten and Evolution layers we have decreased the MAE by approximately 9000. Play around with more layers to see if it's possible to obtain better results.

plot_error(model, X_train, y_train)

Comparing against the validation set

model.evaluate(X_valid, y_valid)
Out: 
{'Explained Variance': 0.8436,
 'MAE': 4861.969,
 'MAPE': 0.1796,
 'MSE': 46412730.8316,
 'RMSE': 6812.689,
 'RMSLE': nan,
 'R2 Score': 0.8047}
plot_error(model, X_valid, y_valid)

model.explain()

Explaining the variance in the Error Plot​

Prior to examining the detailed error plot, it is essential to consider the real-world operational differences among various bulldozer models, as well as the insights provided by subject matter experts (SMEs). These differences are likely to manifest as distinct groupings in the predicted versus actual results. Each model type's unique characteristics—such as age, usage and maintenance history factors that could create these groups, affecting the sale prices and thus the prediction accuracy. Recognizing these potential variances will prepare us to understand and address the disparities in the predictive performance across different Model IDs that the following plot will reveal.

plot_error(model, X_train, y_train, alpha=0.4, color_column="ModelID")

Insights from Scatter Plot Analysis​

The scatter plot displayed above demonstrates a significant variation in the predictive accuracy across different Model IDs, as indicated by the spread of points in relation to the black dashed line, which represents perfect prediction. Models such as those in the yellow cluster are closely aligned with the line, suggesting higher prediction accuracy for these Model IDs. This observation underscores the importance of partitioning the dataset to develop model-specific predictive algorithms. By doing so, we can account for the unique characteristics of each model, which may include factors specific to the model that affect the score contributions.

Creating a Partitioned Model​

from xplainable.core.models import PartitionedRegressor
import pandas as pd
from sklearn.model_selection import train_test_split

# Train your model (this will open an embedded gui)
partitioned_model = PartitionedRegressor(partition_on='ModelID')

# Iterate over the unique values in the partition column
for partition in df.ModelID.value_counts().index[:10].to_list():
    # Get the data for the partition
    part = data[data['ModelID'] == partition].drop(columns=drop_cols)
    x_train_partition, y_train_partition = part.drop('SalePrice', axis=1), part['SalePrice']
    
    # Fit the embedded model
    model_partition = XRegressor()
    model_partition.fit(x_train_partition, y_train_partition, id_columns=id_columns)
    
    #Unhide this if you want to optimise the results
    # network = XEvolutionaryNetwork(model_partition)
    
    # # Add the layers
    # # Start with an initial Tighten layer
    # network.add_layer(
    #     Tighten(
    #         iterations=100,
    #         learning_rate=0.1,
    #         early_stopping=20
    #         )
    #     )
    
    # # Add an Evolve layer with a high severity
    # network.add_layer(
    #     Evolve(
    #         mutations=100,
    #         generations=50,
    #         max_severity=0.5,
    #         max_leaves=20,
    #         early_stopping=20
    #         )
    #     )
    
    # # Add another Evolve layer with a lower severity and reach
    # network.add_layer(
    #     Evolve(
    #         mutations=100,
    #         generations=50,
    #         max_severity=0.3,
    #         max_leaves=15,
    #         early_stopping=20
    #         )
    #     )
    
    # # Add a final Tighten layer with a low learning rate
    # network.add_layer(
    #     Tighten(
    #         iterations=100,
    #         learning_rate=0.025,
    #         early_stopping=20
    #         )
    #     )
    
    # # Fit the network (before or after adding layers)
    # network.fit(X_train.drop(columns=id_columns), y_train)

    # # Run the network
    # network.optimise()

    # Add the model to the partitioned model
    partitioned_model.add_partition(model_partition, partition)

# Predict on the partitioned model
y_pred = partitioned_model.predict(X_valid)
plot_error(partitioned_model, X_train, y_train, color_column="ModelID")

plot_error(partitioned_model, X_valid, y_valid, color_column="ModelID")

Evaluation of Model Predictions Against Validation Data​

The scatter plot illustrates our model's performance on the validation set, comparing the true values against the predicted values for various bulldozer models. While the trend line shows that our model predictions are generally aligned with the true values, there is an observable underprediction across the data points, as evidenced by the mean absolute error (MAE) of 3599 vs 3212 on the train.

Considerations for Model Refinement:​

  • The impact of the mining boom in Australia in 2012, referenced from the Reserve Bank of Australia's report, suggests an economic context that may influence equipment prices. Incorporating macroeconomic indicators could potentially enhance the model's predictive accuracy.

  • Introducing time series features that capture year-over-year changes could offer a more nuanced understanding of price fluctuations over time, rather than relying solely on 'Age at Sale', which may not fully encapsulate such trends.

These considerations point towards the inclusion of external economic factors and more sophisticated time-based features to improve the model's prediction capabilities. Further analysis and iterative model tuning will be required to reduce the prediction error and align the model outputs more closely with the validation data.

Further Investigation:​

  • An analysis of the trend line derived from time series splits (Age at Sale) could reveal insights into future forecasting capabilities. By extending this trend line, we can project forward forecasts that anticipate equipment prices. This approach could be particularly beneficial for capturing the trajectory of market shifts influenced by macroeconomic trends, such as the mining boom.

Should anyone be interested in contributing to the development of this predictive feature or investigating this further, please feel free to add to the issues on our repository or contact us directly at [email protected].

Access model partitions and plot explanations​

partitioned_model.partitions
Out: 
{'4605': <xplainable.core.ml.regression.XRegressor at 0x2d463a1a0>,
 '3538': <xplainable.core.ml.regression.XRegressor at 0x2d45e0160>,
 '4604': <xplainable.core.ml.regression.XRegressor at 0x2d45e1c60>,
 '3170': <xplainable.core.ml.regression.XRegressor at 0x2d45e3af0>,
 '3362': <xplainable.core.ml.regression.XRegressor at 0x2d45e0400>,
 '3537': <xplainable.core.ml.regression.XRegressor at 0x1719da0e0>,
 '4603': <xplainable.core.ml.regression.XRegressor at 0x2d4532290>,
 '3171': <xplainable.core.ml.regression.XRegressor at 0x2d45328f0>,
 '3357': <xplainable.core.ml.regression.XRegressor at 0x171ad8190>,
 '3178': <xplainable.core.ml.regression.XRegressor at 0x1719ca200>}
partitioned_model.partitions['3538'].explain()